不要相信我说了什么，我所说的 未必 是事实，你也 未必 能理解成事实。

总之，经历过回头看就自知了。

莫听穿林打叶声，何妨吟啸且徐行。竹杖芒鞋轻胜马，谁怕，一蓑烟雨任平生。

料峭春风吹酒醒，微冷，山头斜照却相迎。回首向来萧瑟处，归去，也无风雨也无晴。

苏轼，定风波

目录

• toc
  {: toc}

初试

我不记得为什么会做出这些选择，如果非要找个理由。正如Steve Jobs2005年在斯坦福大学演讲时所说的，I follow my heart.

缘起

因为那本叫《Debian标准教程》的书作者曾是byr linux版版主？

看过一本讲TeX的书的作者提到了byr？

然后我来到了byr linux版……

然后接到母亲的电话说可以考北邮，北邮的光学好……我说好啊……

然后准备选通信原理作为专业课，因为我想跨考通信。

买了本书，看了两个月不知道什么玩意……

论坛上有人说，喜欢编程考信安吧，考804信号与系统简单。复试无所谓，大家都不会。

然后不知道为什么我就做了这辈子可能都无法理解的选择……我好像喜欢的不是编程……

复习

八月初说准备考研吧。买本新东方乱序版单词，结果背了一半不到……然后，没有然后了，忘干净了。

开始准备复习高数，李永乐全书还是从学长那要的去年的。这么断续看了一个月，被虐得不轻。好像题目就没有会做且做对的。最后我只看例题了……

专业课一两个月看完了，于是扔掉了……

买本政治吧，大纲解析真贵。我觉得应该最后一个月再看吧，没想到最后根本就没看……

风中劲草买的较早吧，不去自习每天就看看劲草，在寝室读着玩。进度慢爆，拖拖拉拉，同学都复习两遍三遍了我第一遍都没。

打酱油一样的自习，早上八九点起，十一点回，下午三点多去，五点回，晚上不去……一路游击，占个毛线座。

好不容易刻苦了两天，然后生病了。紧接着一个半月书都没碰，天天自己逛公园和看电影去了。

最后只二十来天，想想数学真题看看吧，全书我就看了看例题。历年真题，基本上没答对几道，基本上不是出错就是根本不会……

最后两周，忽然想起貌似英语还没看……把早就买了的真题黄皮书以一天50-180页的速度看完……然后花大概七天背了背附送的单词书，不过貌似没啥用。

考试

运气很好分到个有空调教室，还不靠门。并不是所有人运气都这么好，没空调分阴面的同学在无锡湿冷湿冷的冬天冻爆了。

考前晚上还忙着为恶搞百度百科P图。一副三心二意的模样，考前一周竟然还忙着折腾scipy……

考政治那天晚上漫天大雪，想到要考完试了高兴异常。我又到公园转了一圈才去自习，没怎么看就回寝室，读了十分之一十八大报告就跪了，胡总真不容易。

晚上睡觉时楼上吵得要死，很烦的想带上耳塞等会去掉，结果一戴上睡着了……

第二天下雪了，政治做得没什么特殊感觉，有点冷。

英语，基本都不会，一路瞎蒙，作文划的一塌糊涂。

第二天兴奋的不知两点三点才睡着，想着要结束了就高兴。

上午数学做得史无前例的顺，做完了看还有半个多小时，又翻翻看看，忽然发现有误解的地方，胡乱一改然后就收卷了。然后听旁边同学说数学难……

专业课前半小时后半小时净听监考老师blabla个没完没了，然后简单的会做难得一点看不懂。感觉有的题还是错的。

亢奋的情绪持续了第二天一天，晚上又是两三点才睡着。高兴过头了。

复试完第一件事竟然是装gentoo……太TM屌丝了……

复试

十一号报到十号坐高铁滚到北京，可是奢侈了一把。承蒙同学收留，住在了北科。复试这几天爽爆了，日行二三十公里……北科距北邮40分钟步程(我的走路速度就和慢跑差不多了……)一天至少也得往来一趟。在北邮瞎转一圈，好在北邮小到能从南门看到北门西门看到东门，没觉得太累。第二天又去中科院那边找同学，一路问一路走饶了一大圈又爽了一把，晚上又和同学逛到九点多……但谢谢他(她)们肯花时间招待我，陪我瞎转，很温暖的北京之旅啊。明天打算去北大和清华混吃求陪去，然后哈哈，亦庄。

体检

早上sb了动作慢了一步，同寝的国防生们早冲出去跑操去了。结果面对北科大宿舍新装的伟大门禁我就sb了……等到一会有收垃圾的开了个侧门我终于溜出宿舍了……

结果就是，本打算6点多去校医院排队体检，结果七点才离开北科，又花了四十分钟到北邮校医院。倒是远远看到一长对一路排到教三。然后竟然好多排队不拿表……

拿表，填好，听说可以先体检后缴费就四处找人少的地方。结果后来就sb了，差点没测完……

建议是早去先胸透视力然后抽血，男外别去太晚，其它好像都特快。不然就没人给您体检了……说到男外科我要扯扯，人说北邮校医院有一大绝技——切包皮。我感觉男外那男医生净督促大家切包皮了。

笔试

唯一印象深刻的是我爸发短信说去年要带2b铅笔，让我带上。我根本就不知道，差点悲剧掉。

据说信安的笔试是浮云，我也无心准备。看了看计网，两天内忘光。看看c语言，忘得更快，好歹毕设还和C有些关系。准备看概率论，最后一个字没看……

笔试一个小时，五选三，90道选择题，都是涂答题卡。只要熟悉都不难。

官方说明复试挺重要，由导师自主划线，决定复试人员。但没看到有任何人因为笔试挂掉。

反正没啥好担心的，把复习的时间省下来做其它有意义的事吧。

唯一突然的是，笔试结束后留下来填张表，关于考研成绩，本科专业，是否服从分配，认同团队文化，最早什么时候能来，是否同意调剂到专硕等问题。钮老师还提醒了下复试要带上缴费后打印的报名表，成绩单和简历和刚发的复试表。

面试

面试下午一点半开始，我成绩单简历报名表什么都没打，还放在北科。于是花了四十分钟又从北邮回到北科。半路接到辅导员电话，说查不到入党时间，要问之前负责的C老师，老师听说我在帝都复试，很和善的祝愿我了，瞬间感觉很温暖……回到北科进不去寝室，躺在操场的长椅上看国防生训练，然后就睡着了……

醒来想起入党时间还不知道，打电话问C老师，结果基本是被骂了一遍，呼呼拉拉又打了几通电话，最后老师说，好吧，下午四点半之后打给我吧。结果下午复试时收到老师短信告诉我是几号，真假不知，不过觉得不会坑我，我也不在乎，总之，谢谢C老师。

摸摸兜里就六七块钱了，在北科随便吃了顿饭。然后终于能回寝室了，打印好东西直接又花了四十分钟回到北邮，好在gxy之前带我在北邮逛了一圈，没问几个人就找到新科研楼12层了。

到上面发现我是第二个……qq群里问心是第一个，他进去我就在外面跟一起等待的人扯淡开玩笑。然后我就进去了。

不是什么群面，单独面试感觉挺好。对面一溜五六个老师，我对面坐下。

挺囧的面试……开始让一分钟英文自我介绍，题材随意。我讲了几句看没人听。于是说老师不会的单词用中文行么。老师说算了，你别说了……

然后就开始问我的特长，对什么感兴趣，哪方面研究深入。看我简历说熟悉linux，问我什么是驱动？我完全没明白老师想问什么，很坦然说不知道。出去才想到好像是问我驱动做了什么……也许老师感觉他问了个linux很基本的问题，我都不会，一定在瞎忽悠吧。

发的文章怎么回事，做了多少事，是不是第一作者。项目中你做了多少做了什么……

兴趣点太多，简历写得跳跃度太大。被直接看出来我是很冲动型的人了。不过说我什么坚持不太久我就不高兴了，单linux用这么久了都（好吧，快两年而已）。

问问我成绩怎样，班级排名怎样。班里考研的情况怎样。

问了问我用没用过gcc，用linux能不用gcc？又问用没用过gdb，怎么用，我随口说用过，是的，怎么用gdb调试python我知道，但从来没用gdb调试过c,我说在python中用过(翻译scipy lecture notes时学的)……但老师肯定不知道我的意思……总之，交流能力当时sb了。

然后扯到毕设，好像我扯过去的……扯扯貌似大家没怎么明白，就没怎么问。

问了问为什么选信安，好吧，我就是选信安而已，没啥理由。面对这神回答，老师很无语，我也很无语……

问了我觉得自己有啥优缺点。我说那动手能力和学习能力算有点吧，缺点，有点胆小吧。老师说我面了这么多学生没感觉你胆小，我说我心理有准备。又问问有没其它缺点，我想想说脾气不太好吧，老师说和同学关系不好？我说没有啊，只是有时候感觉太关注自己忽视别人了……

又问如果碰到分到不喜欢的工作会不会很别扭，你这看上去太冲动型了。我说老师您说得对，但我也没啥不感兴趣的，两年半总坚持的了。又说毕竟这个年龄想多看看多学学多折腾折腾，老师点点头，也不知道啥意思。

问我有啥不顺心的事吧，我实在想不到学校有啥，结果就问家里，结果我说了什么来着最后老师竟然来句对不起，我妹的我sb的回了句没事……

问了我到底对啥比较感兴趣，我说数学和计算机吧，渗透性比较强，我也喜欢各种各样交叉的东西，比如说信安。这时候说也有可能让你写文档，我说，好啊，老师您看简历，我很喜欢写东西啊，像我用的archlinux的社区维基上我就写了和翻译了很多东西啊。

被诘难是否没做出什么东西，我说老师您看简历上，做了很多小玩意，虽然不是啥特别好的东西。好歹是我一点一点查资料尝试错误，一点一点实现的。

然后走得时候老师们都很和善，然后我很sb的溜走了，出门被一群人围上了，怎么你这么久，好吧，我感觉就一会的事。

然后……想起来包忘里头了，手里还拿着不知道谁的笔……得，一边跟别人闲扯复试很简单，淡定加油或者开玩笑，一边等下一位出来。下一个出来时很猥琐地进去把包拿走笔放下了……我不知道老师怎么想……这孩子太sb了……

好了，这就是我的研究生复试面试了。sb的人生不需要解释,猥琐的举止无需掩饰……

政面

一早赶去没什么人，老师看看表填好没，随便问问个人信息性格政治倾向一类的。很快。

总结

• 首先，做什么都要有信心吧。不要总是消极自我暗示，要积极自我欺骗。
• 多活动活动，劳逸结合，效率比时间重要。时间是像手中的沙子，抓得越紧，从指缝间流失的越快。
• 考研没什么重要的，身康体健，亲友健在，才是最重要的。

• 别跟大师比较，人比人，比死人。仅仅如一首禅诗：

  To follow the path:（沿着这样一条道路：） look to the master,（寻找大师，） follow the master,（跟随大师，） walk with the master,（与大师通行，） see through the master,（洞察大师，） become the master.（成为大师。）

• 尽人事，听天命。不要胆怯，如果你真的在乎。

• 拿得起，放的下。敢赌敢输，天下之大，何处无容身之地。
• 福兮祸之所倚,祸兮福之所伏。

致谢

• 我的亲人：我不在乎上不上研，只在乎你们。
• 北邮人考研专版：很棒的地方。
• 同学及朋友：谢谢鼓励，谢谢复试时的款待、收留和帮助。
• 老师：从未曾为难过我

总之，很高兴，某件事又结束了。

二维相位去包裹：路径积分法

本文介绍几种基于路径积分的二维相位展开方法。这只是个预告，表示这几天没有闲着。

##为什么二维相位展开会发生问题？

Itoh’s算法的二维扩展

Itoh’s一维相位展开算法很容易扩展到二维。

我们可以先对第一列进行解包裹，然后在此基础上对每行进行解包裹，整个相位图像就都能被解包裹了。

我把扩展的itoh’s算法函数写在util中并导入

1
2

from util import *
from mpl_toolkits.mplot3d import axes3d

获取实验相位函数。我随便选的，应该没问题吧……

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

X, Y, Z = axes3d.get_test_data(0.05)
Z_origin = Z.copy()
plt.figsize(10,10)

plt.subplot(121)
plt.title("origin phase")
plt.imshow(Z, cmap='gray')
ax = plt.gca()
ax.set_axis_off()
colorbar(shrink=.40)

plt.subplot(122)
plt.title("wraped phase")
ax = plt.gca()
ax.set_axis_off()
Z = wrap(Z)
plt.imshow(Z, cmap='gray')
colorbar(shrink=.40)

紧接着用itoh的方法解包裹

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

Z = raster_unwrap2(Z)

plt.figsize(10,10)

plt.subplot(121)
ax = plt.gca()
ax.set_axis_off()
plt.title("unwraped phase")
plt.imshow(Z, cmap='gray')
colorbar(shrink=.40)

plt.subplot(122)
ax = plt.gca()
ax.set_axis_off()
plt.title("erro map")
plt.imshow(Z - Z_origin, cmap='gray')
colorbar(shrink=.40, ticks=None)
np.all(np.abs((Z - Z_origin) < 1e-8))

正如我们所见，这算法效果很好。但是如果数据点出了点问题的话……

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

Z = wrap(Z)

x = np.random.randint(0,120,100)
y = np.random.randint(0,120,100)
Z[x,y] = random.random() * 3

plt.subplot(131)
ax = plt.gca()
ax.set_axis_off()
plt.title("noisy wraped phase")
plt.imshow(Z, cmap='gray')
colorbar(shrink=.25)

Z = raster_unwrap2(Z)

plt.subplot(132)
plt.title("noisy unwraped phase")
ax = plt.gca()
ax.set_axis_off()
plt.imshow(Z, cmap='gray')
colorbar(shrink=.25)

np.all((Z == Z_origin))

plt.subplot(133)
plt.imshow(Z-Z_origin, cmap='gray')
ax = plt.gca()
ax.set_axis_off()
plt.title("error map")
colorbar(shrink=.25)

我加了一些随机噪声，结果就悲剧了。查看误差图给我们了一些启示。在某些出问题的点同一行又边全都出了问题。这种解包裹误差不断扩大并形成一条线的现象，似乎 叫做 拉线 现象(我瞎说的，只是有文献提到过拉线现象，也可能就是这个吧……)。

如果我们对Itoh算法换一种方式扩展。先对第一行解包裹，再在此基础上对各个列解包裹。如果没有有问题的点。结果和第一种扩展方式一样。这给我们个启示，我们可以在任意方向进行相位去包裹，但只要碰到出问题的点都会在之后的积分路径中出现问题。

如果要解包裹出正确的相位，必须绕过这些奇点(有问题的点)。

事实上，相位解包裹领域借用了高等数学和复变函数分析里的理论和观点，并最终形成了相位解包裹问题的留数定理。

相位去包裹留数定理

要想解包裹与路径无关，必须有：

$$\oint \nabla \varphi(r) \cdot dr = 2\pi \times \text{(sum of enclosed residue charges)}$$

在数字图像中，最小的环路是2X2大小的四个像素，如下图所示。通过计算$\sum_{i=1}^4\Delta_i$是否为零，为零则无留数，否则此区域内有正或负留数。实际上，这些有问题的点留数通常不会为0。

为了在任意方向上能进行正确的解包裹，必须想办法使环路内的留数和为0,让积分路径(也是解包裹的路径)绕过这些留数不为0的奇点。

之后所提到的，我所看到的一些路径积分方法。都是通过尽量在可靠的地方进行相位展开(即积分),想办法绕开可能造成问题的点，或平衡这些造成问题的点。

之后做好准备看看那些大牛们发展出的各种路径积分算法吧。我已经佩服的五体投地了。这复杂程度已经让个人Keep it simple， stupid的人生哲学彻底崩坏。各路人马的算法涉及大量数学理论、信号处理理论(mask、noise filter、quality map)、图形形态学操作(开闭操作、连通性判断)、图论(分支切割算法)、最小生成树……完全看跪了……如果你准备动手实践，为了提高时间和空间效率，还有涉及可能是算法设计的东西吧[^1]。

几种路径积分算法

Goldstein 分支切割(branch cut)算法

该算法旨在用最短的路径平衡正负留数，设置解包裹时的障碍，使解包裹路径能按照正确的方向进行。总体上说，该算法相当高效，很容易查看效果。但跟质量图没关系，很容易错误的放置障碍，导致错误的积分路径。

质量引导(quality-guided)算法

先积分质量高的地方，后积分质量低的地方，这基本能保证解包裹路径的正确。这个算法和留数没关系，所以不能保证不出错。但如果有个非常好的质量图，结果相当棒。另外，效率也不错。

掩码分支(mask)算法

前两种算法的杂交，通过某种方法沿着质量最低的地方扩展mask，最后解包裹时绕过这些mask。同上方法，没有好的质量图就别去试。

flynn最小断点算法

很直观的算法。就是把不连续的部分通过增加整数倍个$2\pi$变得连续，虽然说起来没这么简单。该算法用了个生成树，效果很好，还可以用质量图获得更好的效果……

Off-Topic

相位去包裹笔记放在这里。

Footnotes

[^1]:在 Data Analysis with Open Source Tools 一书中推荐过 一本书叫如何设计算法 ，应该值得一看。

一维相位去包裹：原理与仿真

Off-Topic

未觉池塘春草色，阶前梧叶已秋声。

——朱熹

时光飞逝，转眼就要毕业了。前天导师打电话摧我，这几个月就毕业了，你那毕业设计得赶紧做了。想想真是，四年倏忽过去，浑然不觉，已近毕业。

也许这是最后自在平静的日子了，时间却像沙子，越用力抓紧，从指间溜走地越快。

两年过得太快，I follow my passion，做了太多没用的事，有得有失。只是希望：

You can’t connect the dots looking forward you can only connect them looking backwards. So you have to trust that the dots will somehow connect in your future. You have to trust in something: your gut, destiny, life, karma, whatever. Because believing that the dots will connect down the road will give you the confidence to follow your heart, even when it leads you off the well worn path.

——Steve Jobs, Stanford Commencement Adress, 2005

为何相位重要

我大学的专业是光学，时至今日，依然清晰记得两年前在阴暗的实验室做全息摄像的实验。我们几个同学花了好久把光路摆好，让激光打在物体上反光到干版上曝光。然后拿到什么试剂中定形。最后在拿出来用激光将物体的影响重新清晰的放出三维的像时，心里相当兴奋。

为什么全息照相能产生立体的像呢？普通的照相技术都是仅仅记录光波的强度，而不记录相位，因此失去了很多相位中的信息。但全息照相通过相干光(激光)之间的干涉在干版上同时记录下了强度和相位，再用相干光照射干版重放，好像光是从真的物体发出的一样。

如果您也做过这种实验，您应该知道为什么相位重要了。

这还有另一个例子，关于爱因斯坦和蒙娜丽莎。

我从《Two Dimensional Phase Unwrapping Theory Algorithms and Software》中看到了这个例子，然后自己动手用python试了试。

以下代码用来交换两个图片的相位

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>

# <codecell>

from scipy.fftpack import *

# <codecell>

# read file
im_1 = plt.imread('einstein.jpg')
im_2 = plt.imread('monalisa.jpg')

# fft and reverse two images' phase
m_1, p_1 = np.abs(fft2(im_1)), np.angle(fft2(im_1))
m_2, p_2 = np.abs(fft2(im_2)), np.angle(fft2(im_2))

# <codecell>

im_swapphase_1 = np.real(ifft2(m_1 * np.cos(p_2) + m_1 * np.sin(p_2) * 1j))
im_swapphase_2 = np.real(ifft2(m_2 * np.cos(p_1) + m_2 * np.sin(p_1) * 1j))

# <codecell>

plt.figsize(10,10)
plt.gray()
plt.subplot(2,2,1)
plt.imshow(im_1, origin='lower')
plt.subplot(2,2,2)
plt.imshow(im_2, origin='lower')
plt.subplot(2,2,3)
plt.imshow(im_swapphase_1, origin='lower')
plt.subplot(2,2,4)
plt.imshow(im_swapphase_2, origin='lower')

上面两幅图是交换相位前的图，下面两幅是之后的。显然，一团糟，相位中是有信息的。

为何要相位去包裹

简单地说，就是说，任何仪器，比如说量角器，顶多只能测得($-\pi, \pi$]之间的量,但真正的相位角度不该这样，而是分布在实数空间内，应该是测得的($-\pi, \pi$]之间的值的2$\pi$整数倍。

真正的相位值被“包裹”起来了，但为什么要解包裹呢？

不解包裹无法对相位进行计算。

Itoh的路径积分法

我们先定义一个获取包裹相位的算子$\mathscr{W}$，该算子将相位包裹，获取位于$(-\pi,\pi]$之间的包裹相位。

$$\mathscr{W}\varphi = \arctan[\cos(Real \varphi / Img \varphi)]$$

还可以这样写

$$\mathscr{W}{\varphi(n)} = \psi(n) = \varphi(n) + 2\pi k(n)$$

其中$k(n)$是使包裹相位位于$-\pi,\pi$之间的值。

显然包裹相位$\psi(n)$有：

$$\pi \geq \psi(n) \gt -\pi$$

定义差分算子$\Delta$：

$$\Delta {\varphi(n)} = \varphi(n+1) - \varphi(n)$$

$$\Delta {k(n)} = k(n+1) - k(n)$$

计算被包裹相位的差分：

$$\Delta {\psi(n)} = \Delta {\varphi(n)} + 2\pi\Delta {k_1(n)}$$

我们再用$\mathscr{W}$作用于该差分得：

$$\mathscr{W}{\Delta {\psi(n)}} = \Delta {\varphi(n)} + 2\pi[\Delta {k_1(n)} + k_2(n)]$$

显然上式结果应该位于$(-\pi,\pi]$，假如此时还有$\Delta {\varphi(n)}$也位于$(-\pi,\pi]$,则上式右边第二项$2\pi[\Delta {k_1(n) } + k_2(n)]$应该为零，则有：

$$\Delta {\varphi(n)} = \mathscr{W}{\Delta {\psi(n)}} $$

显然，由该差分式可得：

$$\varphi(m) = \varphi(0) + \sum_{n=0}^{m-1} \mathscr{W}{\Delta {\mathscr{W}{\varphi(n)}}} $$

上式说明，真实相位可以通过对包裹相位的差分的包裹进行积分求得。

于是itoh的一维相位去包裹算法综述如下：

对信号相位数组$\psi(i),0 \leq i \leq N-1$

• 计算相位差分$D(i) = \psi(i+1)-\psi(i), i=0,\ldots,N-2$
• 计算包裹的相位差分$\Delta(i) = \Delta{D(i)}, i=0,\ldots,N-2$
• 初始化初值$\varphi(0)= \psi(0)$
• 累加解包裹$\varphi(i) = \psi(i) + \Delta(i)$

Itoh的方法很简单实用，但受到两个重要因素的影响：相位失真和噪声。下面仿真两种情况的影响。

仿真

对正弦波相位函数(间谐波，一切波的基础)

$$\varphi(t) = 10\sin(10t), 0 \leq t \leq 1$$

通过计算可以得知，使之不产生相位失真，区间内至少有32个采样点:

对相位变化有

$$\Delta \varphi = \dot{\varphi}\Delta t$$

其中 $\dot{\varphi} = d\varphi / dt = 100\cos10t$ ，可看出相位在 $n\pi/10$ 取得极值，加上条件：

$$\left\vert\Delta \varphi\right\vert \lt \pi$$

则

$$\left\vert \frac{100}{N} \cos 10t \right\vert \lt \pi$$

得

$$N \gt 31.83$$

则对这个函数如果采样率低于32就会产生相位失真。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>

# <codecell>

import numpy as np
import matplotlib.pyplot as plt

# <codecell>

# phase function
def pf(x):
    return 10 * np.sin(10 * x)
# wrap function
def wrap(x):
    return np.arctan2(np.sin(x), np.cos(x))
# unwrap function
def wrap_diff(x):
     return wrap(np.diff(x))

def unwrap(x):
    y = x
    y[0] = x[0]
    for i in range(len(x) - 1):
        i += 1
        y[i] = y[i - 1] + wrap_diff(x)[i - 1]
    return np.array(y)

def noise(x, snr):
    return x + np.random.normal(loc=0.0, scale = np.sqrt(np.max(x) / 2 ** snr), size=len(x))

def show_unwrap(x,y,y_w,t,p):
    plt.ylim((-12,12))
    plt.ylabel("Phase in Radians")
    p_w, = plt.plot(x, y_w,'o:')
    p_o, = plt.plot(t,p,':')
    p_u, = plt.plot(x,y,'s')
    plt.legend([p_w, p_o, p_u], ["sampled wrapped phase", "original phase function", "unwrapped phase"])
    

# <headingcell level=1>

# 当取样点为num时

# <codecell>

# Origin
t = np.arange(0,1,0.01)
p = pf(t)

# plot
## num 50
plt.subplot(311)
plt.title("num = 50")
# sampled data
num = 50
x = np.linspace(0,1,num)
# wrapped phase
y_w = wrap(pf(x))
# unwrapped wraped phase
y = unwrap(wrap(pf(x)))
show_unwrap(x, y, y_w, t, p)
## num = 32
plt.subplot(312)
plt.title("num = 32")
# sampled data
num = 32
x = np.linspace(0,1,num)
# wrapped phase
y_w = wrap(pf(x))
# unwrapped wraped phase
y = unwrap(wrap(pf(x)))
show_unwrap(x, y, y_w, t, p)
## num = 31
plt.subplot(313)
plt.title("num = 31")
# sampled data
num = 31
x = np.linspace(0,1,num)
# wrapped phase
y_w = wrap(pf(x))
# unwrapped wraped phase
y = unwrap(wrap(pf(x)))
show_unwrap(x, y, y_w, t, p)
plt.xlabel("Relative Time")


# <headingcell level=1>

# 噪声

# <codecell>

# noise influence
x = np.linspace(0,1,200)
y = pf(x)
y_10 = unwrap(wrap(noise(y,10))) + 20
y_5 = unwrap(noise(wrap(y), 5)) + 40
y_2 = unwrap(noise(wrap(y), 2)) + 60
y_1 = unwrap(noise(wrap(y), 1)) + 80
# plot
plt.xlabel("Relative Time")
plt.ylabel("Phase in Radians")
p_o, = plt.plot(t,p)
p_1, = plt.plot(x,y_1,'--')
p_2, = plt.plot(x,y_2,'-.')
p_5, = plt.plot(x,y_5,':')
p_10, = plt.plot(x,y_10,':')
plt.legend([p_o, p_1, p_2, p_5, p_10], ["Origin", "SNR=1", "SNR=2", "SNR=5", "SNR=10"])

# <codecell>

我们可以看到：

• 在采样率低于某一关键值32时基本没法解出正确相位。

  

• 信噪比越低，解包裹效果越差。
  

Summary

如果是二维相位去包裹问题，还有个奇点的问题。噪音、相位失真、奇点，似乎是所有相位解缠算法必须面对的三大问题。

相关Ipython Notebook文件和Einstein和Monalisa的图像在此下载

高级Python结构

原谅渣翻译，可能仅仅是给我自己看的。本来多年(也就几个月吧)之前将此文投递到OSChina翻译频道结果被以代码太多为由被拒，于是译者只好用自己的渣英语渣水平翻译给自己看了……，期待各路大婶们指正……

翻译自http://scipy-lectures.github.com/advanced/advanced_python/index.html

作者: Zbigniew Jędrzejewski-Szmek

这章有关Python中被认为高级的特性——就是说并不是每个语言都有的，也是说它们可能在更复杂的程序或库中更有用，但不是说特别特殊或特别复杂。

强调这点很重要：这一章仅仅关于语言自身——关于辅之以Python的标准库功能的特殊语法所支持的特性，不包括那些智能的外部模块实现。

在开发Python程序语言的过程中，它的语法，独一无二。因为它非常透明。建议的更改通过不同的角度评估并在公开邮件列表讨论，最终决定考虑到假设用例的重要性、添加更多特性的负担，其余语法的一致性、是否建议的变种易于读写和理解之间的平衡。这个过程由Python Enhancement Proposals(PEPs)的形式规范。最终这一章节中描述的特性在证明它们确实解决实际问题并且使用起来尽可能简单后被添加。

目录

• toc
  {: toc}

迭代器(Iterators), 生成表达式(generator expressions)和生成器(generators)

迭代器

简单

重复工作是浪费，将不同“土生土长”的方法替换为标准特性换来的是更加易于阅读和操作。

Guido van Rossum — Adding Optional Static Typing to Python

迭代器是依附于迭代协议的对象——基本意味它有一个next方法(method)，当调用时，返回序列中的下一个项目。当无项目可返回时，引发(raise)StopIteration异常。

迭代对象允许一次循环。它保留单次迭代的状态(位置)，或从另一个角度讲，每次循环序列都需要一个迭代对象。这意味我们可以同时迭代同一个序列不只一次。将迭代逻辑和序列分离使我们有更多的迭代方式。

调用一个容器(container)的__iter__方法创建迭代对象是掌握迭代器最直接的方式。iter函数为我们节约一些按键。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

>>> nums = [1,2,3]      # note that ... varies: these are different objects
>>> iter(nums)                           
<listiterator object at ...>
>>> nums.__iter__()                      
<listiterator object at ...>
>>> nums.__reversed__()                  
<listreverseiterator object at ...>

>>> it = iter(nums)
>>> next(it)            # next(obj) simply calls obj.next()
1
>>> it.next()
2
>>> next(it)
3
>>> next(it)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
StopIteration

当在循环中使用时，StopIteration被接受并停止循环。但通过显式引发(invocation)，我们看到一旦迭代器元素被耗尽，存取它将引发异常。

使用for...in循环也使用__iter__方法。这允许我们透明地开始对一个序列迭代。但是如果我们已经有一个迭代器，我们想在for循环中能同样地使用它们。为了实现这点，迭代器除了next还有一个方法__iter__来返回迭代器自身(self)。

Python中对迭代器的支持无处不在：标准库中的所有序列和无序容器都支持。这个概念也被拓展到其它东西：例如file对象支持行的迭代。

1
2
3

>>> f = open('/etc/fstab')
>>> f is f.__iter__()
True

file自身就是迭代器，它的__iter__方法并不创建一个单独的对象：仅仅单线程的顺序读取被允许。

生成表达式

第二种创建迭代对象的方式是通过 生成表达式(generator expression) ，列表推导(list comprehension)的基础。为了增加清晰度，生成表达式总是封装在括号或表达式中。如果使用圆括号，则创建了一个生成迭代器(generator iterator)。如果是方括号，这一过程被‘短路’我们获得一个列表list。

1
2
3
4
5
6

>>> (i for i in nums)                    
<generator object <genexpr> at 0x...>
>>> [i for i in nums]
[1, 2, 3]
>>> list(i for i in nums)
[1, 2, 3]

在Python 2.7和 3.x中列表表达式语法被扩展到 字典和集合表达式。一个集合set当生成表达式是被大括号封装时被创建。一个字典dict在表达式包含key:value形式的键值对时被创建：

1
2
3
4

>>> {i for i in range(3)}   
set([0, 1, 2])
>>> {i:i**2 for i in range(3)}   
{0: 0, 1: 1, 2: 4}

如果您不幸身陷古老的Python版本中，这个语法有点糟：

1
2
3
4

>>> set(i for i in 'abc')
set(['a', 'c', 'b'])
>>> dict((i, ord(i)) for i in 'abc')
{'a': 97, 'c': 99, 'b': 98}

生成表达式相当简单，不用多说。只有一个陷阱值得提及：在版本小于3的Python中索引变量(i)会泄漏。

生成器

生成器

生成器是产生一列结果而不是单一值的函数。

David Beazley — A Curious Course on Coroutines and Concurrency

第三种创建迭代对象的方式是调用生成器函数。一个 生成器(generator) 是包含关键字yield的函数。值得注意，仅仅是这个关键字的出现完全改变了函数的本质：yield语句不必引发(invoke)，甚至不必可接触。但让函数变成了生成器。当一个函数被调用时，其中的指令被执行。而当一个生成器被调用时，执行在其中第一条指令之前停止。生成器的调用创建依附于迭代协议的生成器对象。就像常规函数一样，允许并发和递归调用。

当next被调用时，函数执行到第一个yield。每次遇到yield语句获得一个作为next返回的值，在yield语句执行后，函数的执行又被停止。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

>>> def f():
...   yield 1
...   yield 2
>>> f()                                   
<generator object f at 0x...>
>>> gen = f()
>>> gen.next()
1
>>> gen.next()
2
>>> gen.next()
Traceback (most recent call last):
 File "<stdin>", line 1, in <module>
StopIteration

让我们遍历单个生成器函数调用的整个历程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

>>> def f():
...   print("-- start --")
...   yield 3
...   print("-- middle --")
...   yield 4
...   print("-- finished --")
>>> gen = f()
>>> next(gen)
-- start --
3
>>> next(gen)
-- middle --
4
>>> next(gen)                            
-- finished --
Traceback (most recent call last):
 ...
StopIteration

相比常规函数中执行f()立即让print执行，gen不执行任何函数体中语句就被赋值。只有当gen.next()被next调用，直到第一个yield部分的语句才被执行。第二个语句打印-- middle --并在遇到第二个yield时停止执行。第三个next打印-- finished --并且到函数末尾，因为没有yield，引发了异常。

当函数yield之后控制返回给调用者后发生了什么？每个生成器的状态被存储在生成器对象中。从这点看生成器函数，好像它是运行在单独的线程，但这仅仅是假象：执行是严格单线程的，但解释器保留和存储在下一个值请求之间的状态。^1

为何生成器有用？正如关于迭代器这部分强调的，生成器函数只是创建迭代对象的又一种方式。一切能被yield语句完成的东西也能被next方法完成。然而，使用函数让解释器魔力般地创建迭代器有优势。一个函数可以比需要next和__iter__方法的类定义短很多。更重要的是，相比不得不对迭代对象在连续next调用之间传递的实例(instance)属性来说，生成器的作者能更简单的理解局限在局部变量中的语句。

还有问题是为何迭代器有用？当一个迭代器用来驱动循环，循环变得简单。迭代器代码初始化状态，决定是否循环结束，并且找到下一个被提取到不同地方的值。这凸显了循环体——最值得关注的部分。除此之外，可以在其它地方重用迭代器代码。

双向通信

每个yield语句将一个值传递给调用者。这就是为何PEP 255引入生成器(在Python2.2中实现)。但是相反方向的通信也很有用。一个明显的方式是一些外部(extern)语句，或者全局变量或共享可变对象。通过将先前无聊的yield语句变成表达式，直接通信因PEP 342成为现实(在2.5中实现)。当生成器在yield语句之后恢复执行时，调用者可以对生成器对象调用一个方法，或者传递一个值 给 生成器，然后通过yield语句返回，或者通过一个不同的方法向生成器注入异常。

第一个新方法是send(value)，类似于next()，但是将value传递进作为yield表达式值的生成器中。事实上，g.next()和g.send(None)是等效的。

第二个新方法是throw(type, value=None, traceback=None)，等效于在yield语句处

1

raise type, value, traceback

不像raise(从执行点立即引发异常)，throw()首先恢复生成器，然后仅仅引发异常。选用单次throw就是因为它意味着把异常放到其它位置，并且在其它语言中与异常有关。

当生成器中的异常被引发时发生什么？它可以或者显式引发，当执行某些语句时可以通过throw()方法注入到yield语句中。任一情况中，异常都以标准方式传播：它可以被except和finally捕获，或者造成生成器的中止并传递给调用者。

因完整性缘故，值得提及生成器迭代器也有close()方法，该方法被用来让本可以提供更多值的生成器立即中止。它用生成器的__del__方法销毁保留生成器状态的对象。

让我们定义一个只打印出通过send和throw方法所传递东西的生成器。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

>>> import itertools
>>> def g():
...     print '--start--'
...     for i in itertools.count():
...         print '--yielding %i--' % i
...         try:
...             ans = yield i
...         except GeneratorExit:
...             print '--closing--'
...             raise
...         except Exception as e:
...             print '--yield raised %r--' % e
...         else:
...             print '--yield returned %s--' % ans

>>> it = g()
>>> next(it)
--start--
--yielding 0--
0
>>> it.send(11)
--yield returned 11--
--yielding 1--
1
>>> it.throw(IndexError)
--yield raised IndexError()--
--yielding 2--
2
>>> it.close()
--closing--

注意： next还是__next__?

在Python 2.x中，接受下一个值的迭代器方法是next，它通过全局函数next显式调用，意即它应该调用__next__。就像全局函数iter调用__iter__。这种不一致在Python 3.x中被修复，it.next变成了it.__next__。对于其它生成器方法——send和throw情况更加复杂，因为它们不被解释器隐式调用。然而，有建议语法扩展让continue带一个将被传递给循环迭代器中send的参数。如果这个扩展被接受，可能gen.send会变成gen.__send__。最后一个生成器方法close显然被不正确的命名了，因为它已经被隐式调用。

链式生成器

注意： 这是PEP 380的预览(还未被实现，但已经被Python3.3接受)[^2]

比如说我们正写一个生成器，我们想要yield一个第二个生成器——一个子生成器(subgenerator)——生成的数。如果仅考虑产生(yield)的值，通过循环可以不费力的完成：

1
2
3

subgen = some_other_generator()
for v in subgen:
    yield v

然而，如果子生成器需要调用send()、throw()和close()和调用者适当交互的情况下，事情就复杂了。yield语句不得不通过类似于前一章节部分定义的try...except...finally结构来保证“调试”生成器函数。这种代码在PEP 380中提供，现在足够拿出将在Python 3.3中引入的新语法了：

1

yield from some_other_generator()

像上面的显式循环调用一样，重复从some_other_generator中产生值直到没有值可以产生，但是仍然向子生成器转发send、throw和close。

装饰器

总结

这个语言中令人激动的特性几乎充满歉意的，考虑到它可能没这么有用。

Bruce Eckel — An Introduction to Python Decorators

因为函数或类都是对象，它们也能被四处传递。它们又是可变对象，可以被更改。在函数或类对象创建后但绑定到名字前更改之的行为为装饰(decorator)。^4

“装饰器”后隐藏了两种意思——一是函数起了装饰作用，例如，执行真正的工作，另一个是依附于装饰器语法的表达式，例如，at符号和装饰函数的名称。

函数可以通过函数装饰器语法装饰：

1
2
3

@decorator             # ②
def function():        # ①
    pass

• 函数以标准方式定义。①
• 以@做为定义为装饰器函数前缀的表达式②。在 @ 后的部分必须是简单的表达式，通常只是函数或类的名字。这一部分先求值，在下面的定义的函数准备好后，装饰器被新定义的函数对象作为单个参数调用。装饰器返回的值附着到被装饰的函数名。

装饰器可以应用到函数和类上。对类语义很明晰——类定义被当作参数来调用装饰器，无论返回什么都赋给被装饰的名字。

在装饰器语法实现前(PEP 318)，通过将函数和类对象赋给临时变量然后显式调用装饰器然后将返回值赋给函数名，可以完成同样的事。这似乎要打更多的字，也确实装饰器函数名用了两次同时临时变量要用至少三次，很容易出错。以上实例相当于：

1
2
3

def function():                  # ①
    pass
function = decorator(function)   # ②

装饰器可以堆栈(stacked)——应用的顺序是从底到上或从里到外。就是说最初的函数被当作第一次参数器的参数，无论返回什么都被作为第二个装饰器的参数……无论最后一个装饰器返回什么都被依附到最初函数的名下。

装饰器语法因其可读性被选择。因为装饰器在函数头部前被指定，显然不是函数体的一部分，它只能对整个函数起作用。以@为前缀的表达式又让它明显到不容忽视(根据PEP叫在您脸上……：）)。当多个装饰器被应用时，每个放在不同的行非常易于阅读。

代替和调整原始对象

装饰器可以或者返回相同的函数或类对象或者返回完全不同的对象。第一种情况中，装饰器利用函数或类对象是可变的添加属性，例如向类添加文档字符串(docstring).装饰器甚至可以在不改变对象的情况下做有用的事，例如在全局注册表中注册装饰的类。在第二种情况中，简直无所不能：当什么不同的东西取代了被装饰的类或函数，新对象可以完全不同。然而这不是装饰器的目的：它们意在改变装饰对象而非做不可预料的事。因此当一个函数在装饰时被完全替代成不同的函数时，新函数通常在一些准备工作后调用原始函数。同样，当一个类被装饰成一个新类时，新类通常源于被装饰类。当装饰器的目的是“每次都”做什么，像记录每次对被装饰函数的调用，只有第二类装饰器可用。另一方面，如果第一类足够了，最好使用它因为更简单。^3

实现类和函数装饰器

对装饰器惟一的要求是它能够单参数调用。这意味着装饰器可以作为常规函数或带有__call__方法的类的实现，理论上，甚至lambda函数也行。

让我们比较函数和类方法。装饰器表达式(@后部分)可以只是名字。只有名字的方法很好(打字少，看起来整洁等)，但是只有当无需用参数定制装饰器时才可能。被写作函数的装饰器可以用以下两种方式：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

>>> def simple_decorator(function):
...   print "doing decoration"
...   return function
>>> @simple_decorator
... def function():
...   print "inside function"
doing decoration
>>> function()
inside function

>>> def decorator_with_arguments(arg):
...   print "defining the decorator"
...   def _decorator(function):
...       # in this inner function, arg is available too
...       print "doing decoration,", arg
...       return function
...   return _decorator
>>> @decorator_with_arguments("abc")
... def function():
...   print "inside function"
defining the decorator
doing decoration, abc
>>> function()
inside function

这两个装饰器属于返回被装饰函数的类别。如果它们想返回新的函数，需要额外的嵌套，最糟的情况下，需要三层嵌套。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

>>> def replacing_decorator_with_args(arg):
...   print "defining the decorator"
...   def _decorator(function):
...       # in this inner function, arg is available too
...       print "doing decoration,", arg
...       def _wrapper(*args, **kwargs):
...           print "inside wrapper,", args, kwargs
...           return function(*args, **kwargs)
...       return _wrapper
...   return _decorator
>>> @replacing_decorator_with_args("abc")
... def function(*args, **kwargs):
...     print "inside function,", args, kwargs
...     return 14
defining the decorator
doing decoration, abc
>>> function(11, 12)
inside wrapper, (11, 12) {}
inside function, (11, 12) {}
14

_wrapper函数被定义为接受所有位置和关键字参数。通常我们不知道哪些参数被装饰函数会接受，所以wrapper将所有东西都创递给被装饰函数。一个不幸的结果就是显式参数很迷惑人。

相比定义为函数的装饰器，定义为类的复杂装饰器更简单。当对象被创建，__init__方法仅仅允许返回None，创建的对象类型不能更改。这意味着当装饰器被定义为类时，使用无参数的形式没什么意义：最终被装饰的对象只是装饰类的一个实例而已，被构建器(constructor)调用返回，并不非常有用。讨论在装饰表达式中给出参数的基于类的装饰器，__init__方法被用来构建装饰器。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

>>> class decorator_class(object):
...   def __init__(self, arg):
...       # this method is called in the decorator expression
...       print "in decorator init,", arg
...       self.arg = arg
...   def __call__(self, function):
...       # this method is called to do the job
...       print "in decorator call,", self.arg
...       return function
>>> deco_instance = decorator_class('foo')
in decorator init, foo
>>> @deco_instance
... def function(*args, **kwargs):
...   print "in function,", args, kwargs
in decorator call, foo
>>> function()
in function, () {}

相对于正常规则(PEP 8)由类写成的装饰器表现得更像函数，因此它们的名字以小写字母开始。

事实上，创建一个仅返回被装饰函数的新类没什么意义。对象应该有状态，这种装饰器在装饰器返回新对象时更有用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

>>> class replacing_decorator_class(object):
...   def __init__(self, arg):
...       # this method is called in the decorator expression
...       print "in decorator init,", arg
...       self.arg = arg
...   def __call__(self, function):
...       # this method is called to do the job
...       print "in decorator call,", self.arg
...       self.function = function
...       return self._wrapper
...   def _wrapper(self, *args, **kwargs):
...       print "in the wrapper,", args, kwargs
...       return self.function(*args, **kwargs)
>>> deco_instance = replacing_decorator_class('foo')
in decorator init, foo
>>> @deco_instance
... def function(*args, **kwargs):
...   print "in function,", args, kwargs
in decorator call, foo
>>> function(11, 12)
in the wrapper, (11, 12) {}
in function, (11, 12) {}

像这样的装饰器可以做任何事，因为它能改变被装饰函数对象和参数，调用被装饰函数或不调用，最后改变返回值。

复制原始函数的文档字符串和其它属性

当新函数被返回代替装饰前的函数时，不幸的是原函数的函数名，文档字符串和参数列表都丢失了。这些属性可以部分通过设置__doc__(文档字符串)，__module__和__name__(函数的全称)、__annotations__(Python 3中关于参数和返回值的额外信息)移植到新函数上，这些工作可通过functools.update_wrapper自动完成。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

>>> import functools
>>> def better_replacing_decorator_with_args(arg):
...   print "defining the decorator"
...   def _decorator(function):
...       print "doing decoration,", arg
...       def _wrapper(*args, **kwargs):
...           print "inside wrapper,", args, kwargs
...           return function(*args, **kwargs)
...       return functools.update_wrapper(_wrapper, function)
...   return _decorator
>>> @better_replacing_decorator_with_args("abc")
... def function():
...     "extensive documentation"
...     print "inside function"
...     return 14
defining the decorator
doing decoration, abc
>>> function                           
<function function at 0x...>
>>> print function.__doc__
extensive documentation

一件重要的东西是从可迁移属性列表中所缺少的：参数列表。参数的默认值可以通过__defaults__、__kwdefaults__属性更改，但是不幸的是参数列表本身不能被设置为属性。这意味着help(function)将显式无用的参数列表，使使用者迷惑不已。一个解决此问题有效但是丑陋的方式是使用eval动态创建wrapper。可以使用外部external模块自动实现。它提供了对decorator装饰器的支持，该装饰器接受wrapper并将之转换成保留函数签名的装饰器。

综上，装饰器应该总是使用functools.update_wrapper或者其它方式赋值函数属性。

标准库中的示例

首先要提及的是标准库中有一些实用的装饰器，有三种装饰器：

• classmethod让一个方法变成“类方法”，即它能够无需创建实例调用。当一个常规方法被调用时，解释器插入实例对象作为第一个参数self。当类方法被调用时，类本身被给做第一个参数，一般叫cls。

  类方法也能通过类命名空间读取，所以它们不必污染模块命名空间。类方法可用来提供替代的构建器(constructor):

  1 2 3 4 5 6 7 8

  class Array(object): def __init__(self, data): self.data = data @classmethod def fromfile(cls, file): data = numpy.load(file) return cls(data)

这比用一大堆标记的__init__简单多了。

• staticmethod应用到方法上让它们“静态”，例如，本来一个常规函数，但通过类命名空间存取。这在函数仅在类中需要时有用(它的名字应该以_为前缀)，或者当我们想要用户以为方法连接到类时也有用——虽然对实现本身不必要。

• property是对getter和setter问题Python风格的答案。通过property装饰的方法变成在属性存取时自动调用的getter。

  1 2 3 4 5 6 7 8 9

  >>> class A(object): ... @property ... def a(self): ... "an important attribute" ... return "a value" >>> A.a <property object at 0x...> >>> A().a 'a value'

  例如A.a是只读属性，它已经有文档了：help(A)包含从getter方法获取的属性a的文档字符串。将a定义为property使它能够直接被计算，并且产生只读的副作用，因为没有定义任何setter。

  为了得到setter和getter，显然需要两个方法。从Python 2.6开始首选以下语法：

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

  class Rectangle(object): def __init__(self, edge): self.edge = edge @property def area(self): """Computed area. Setting this updates the edge length to the proper value. """ return self.edge**2 @area.setter def area(self, area): self.edge = area ** 0.5

通过`property`装饰器取代带一个属性(property)对象的getter方法，以上代码起作用。这个对象反过来有三个可用于装饰器的方法`getter`、`setter`和`deleter`。它们的作用就是设定属性对象的getter、setter和deleter(被存储为`fget`、`fset`和`fdel`属性(attributes))。当创建对象时，getter可以像上例一样设定。当定义setter时，我们已经在`area`中有property对象，可以通过`setter`方法向它添加setter，一切都在创建类时完成。

之后，当类实例创建后，property对象和特殊。当解释器执行属性存取、赋值或删除时，其执行被下放给property对象的方法。

为了让一切一清二楚[^5]，让我们定义一个“调试”例子：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
>>> class D(object):
...    @property
...    def a(self):
...      print "getting", 1
...      return 1
...    @a.setter
...    def a(self, value):
...      print "setting", value
...    @a.deleter
...    def a(self):
...      print "deleting"
>>> D.a                                    
<property object at 0x...>
>>> D.a.fget                               
<function a at 0x...>
>>> D.a.fset                               
<function a at 0x...>
>>> D.a.fdel                               
<function a at 0x...>
>>> d = D()               # ... varies, this is not the same `a` function
>>> d.a
getting 1
1
>>> d.a = 2
setting 2
>>> del d.a
deleting
>>> d.a
getting 1
1



属性(property)是对装饰器语法的一点扩展。使用装饰器的一大前提——命名不重复——被违反了，但是目前没什么更好的发明。为getter，setter和deleter方法使用相同的名字还是个好的风格。

一些其它更新的例子包括：

• functools.lru_cache记忆任意维持有限 参数：结果 对的缓存函数(Python 3.2)
• functools.total_ordering是一个基于单个比较方法而填充丢失的比较(ordering)方法(__lt__,__gt__，__le__等等)的类装饰器。

函数的废弃

比如说我们想在第一次调用我们不希望被调用的函数时在标准错误打印一个废弃函数警告。如果我们不像更改函数，我们可用装饰器

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

class deprecated(object):
    """Print a deprecation warning once on first use of the function.

    >>> @deprecated()                    # doctest: +SKIP
    ... def f():
    ...     pass
    >>> f()                              # doctest: +SKIP
    f is deprecated
    """
    def __call__(self, func):
        self.func = func
        self.count = 0
        return self._wrapper
    def _wrapper(self, *args, **kwargs):
        self.count += 1
        if self.count == 1:
            print self.func.__name__, 'is deprecated'
        return self.func(*args, **kwargs)

也可以实现成函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

def deprecated(func):
    """Print a deprecation warning once on first use of the function.

    >>> @deprecated                      # doctest: +SKIP
    ... def f():
    ...     pass
    >>> f()                              # doctest: +SKIP
    f is deprecated
    """
    count = [0]
    def wrapper(*args, **kwargs):
        count[0] += 1
        if count[0] == 1:
            print func.__name__, 'is deprecated'
        return func(*args, **kwargs)
    return wrapper

while-loop移除装饰器

例如我们有个返回列表的函数，这个列表由循环创建。如果我们不知道需要多少对象，实现这个的标准方法如下：

1
2
3
4
5
6
7
8

def find_answers():
    answers = []
    while True:
        ans = look_for_next_answer()
        if ans is None:
            break
        answers.append(ans)
    return answers

只要循环体很紧凑，这很好。一旦事情变得更复杂，正如真实的代码中发生的那样，这就很难读懂了。我们可以通过yield语句简化它，但之后用户不得不显式调用嗯list(find_answers())。

我们可以创建一个为我们构建列表的装饰器：

1
2
3
4

def vectorized(generator_func):
    def wrapper(*args, **kwargs):
        return list(generator_func(*args, **kwargs))
    return functools.update_wrapper(wrapper, generator_func)

然后函数变成这样：

1
2
3
4
5
6
7

@vectorized
def find_answers():
    while True:
        ans = look_for_next_answer()
        if ans is None:
            break
        yield ans

插件注册系统

这是一个仅仅把它放进全局注册表中而不更改类的类装饰器，它属于返回被装饰对象的装饰器。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

class WordProcessor(object):
    PLUGINS = []
    def process(self, text):
        for plugin in self.PLUGINS:
            text = plugin().cleanup(text)
        return text

    @classmethod
    def plugin(cls, plugin):
        cls.PLUGINS.append(plugin)

@WordProcessor.plugin
class CleanMdashesExtension(object):
    def cleanup(self, text):
        return text.replace('—', u'\N{em dash}')

这里我们使用装饰器完成插件注册。我们通过一个名词调用装饰器而不是一个动词，因为我们用它来声明我们的类是WordProcessor的一个插件。plugin方法仅仅将类添加进插件列表。

关于插件自身说下：它用真正的Unicode中的破折号符号替代HTML中的破折号。它利用unicode literal notation通过它在unicode数据库中的名称(“EM DASH”)插入一个符号。如果直接插入Unicode符号，将不可能区分所插入的和源程序中的破折号。

更多例子和参考

• PEP 310(函数和方法装饰语法)
• PEP 3129(类装饰语法)
• http://wiki.python.org/moin/PythonDecoratorLibrary
• http://docs.python.org/dev/library/functools.html
• http://pypi.python.org/pypi/decorator
• Bruce Eckel
  • 装饰器I:介绍Python装饰器
  • Python装饰器II：装饰器参数
  • Python装饰器III：一个基于装饰器构建的系统

上下文管理器

上下文管理器是可以在with语句中使用，拥有__enter__和__exit__方法的对象。

1
2

with manager as var:
    do_something(var)

相当于以下情况的简化：

1
2
3
4
5

var = manager.__enter__()
try:
    do_something(var)
finally:
    manager.__exit__()

换言之，PEP 343中定义的上下文管理器协议允许将无聊的try...except...finally结构抽象到一个单独的类中，仅仅留下关注的do_something部分。

1.__enter__方法首先被调用。它可以返回赋给var的值。as部分是可选的：如果它不出现，enter的返回值简单地被忽略。
2.with语句下的代码被执行。就像try子句，它们或者成功执行到底，或者break，continue或return，或者可以抛出异常。无论哪种情况，该块结束后，__exit__方法被调用。如果抛出异常，异常信息被传递给__exit__，这将在下一章节讨论。通常情况下，异常可被忽略，就像在finally子句中一样，并且将在__exit__结束后重新抛出。

比如说我们想确认一个文件在完成写操作之后被立即关闭：

1
2
3
4
5
6
7
8
9

>>> class closing(object):
...   def __init__(self, obj):
...     self.obj = obj
...   def __enter__(self):
...     return self.obj
...   def __exit__(self, *args):
...     self.obj.close()
>>> with closing(open('/tmp/file', 'w')) as f:
...   f.write('the contents\n')

这里我们确保了当with块退出时调用了f.close()。因为关闭文件是非常常见的操作，该支持已经出现在file类之中。它有一个__exit__方法调用close，并且本身可作为上下文管理器。

1
2

>>> with open('/tmp/file', 'a') as f:
...   f.write('more contents\n')

try...finally常见的用法是释放资源。各种不同的情况实现相似：在__enter__阶段资源被获得，在__exit__阶段释放，如果抛出异常也被传递。正如文件操作，往往这是对象使用后的自然操作，内置支持使之很方便。每一个版本，Python都在更多的地方提供支持。

• 所有类似文件的对象：
  • file ➔ 自动关闭
  • fileinput,tempfile(py >= 3.2)
  • bz2.BZ2File，gzip.GzipFile, tarfile.TarFile,zipfile.ZipFile
  • ftplib, nntplib ➔ 关闭连接(py >= 3.2)
• 锁
  • multiprocessing.RLock ➔ 锁定和解锁
  • multiprocessing.Semaphore
  • memoryview ➔ 自动释放(py >= 3.2 或 3.3)
• decimal.localcontext➔ 暂时更改计算精度
• _winreg.PyHKEY ➔ 打开和关闭Hive Key
• warnings.catch_warnings ➔ 暂时杀死(kill)警告
• contextlib.closing ➔ 如上例，调用close
• 并行编程
  • concurrent.futures.ThreadPoolExecutor ➔ 并行调用然后杀掉线程池(py >= 3.2)
  • concurrent.futures.ProcessPoolExecutor ➔ 并行调用并杀死进程池(py >= 3.2)
  • nogil ➔ 暂时解决GIL问题(仅仅cyphon ：（)

捕获异常

当一个异常在with块中抛出时，它作为参数传递给__exit__。三个参数被使用，和sys.exc_info()返回的相同：类型、值和回溯(traceback)。当没有异常抛出时，三个参数都是None。上下文管理器可以通过从__exit__返回一个真(True)值来“吞下”异常。例外可以轻易忽略，因为如果__exit__不使用return直接结束，返回None——一个假(False)值，之后在__exit__结束后重新抛出。

捕获异常的能力创造了有意思的可能性。一个来自单元测试的经典例子——我们想确保一些代码抛出正确种类的异常：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

class assert_raises(object):
    # based on pytest and unittest.TestCase
    def __init__(self, type):
        self.type = type
    def __enter__(self):
        pass
    def __exit__(self, type, value, traceback):
        if type is None:
            raise AssertionError('exception expected')
        if issubclass(type, self.type):
            return True # swallow the expected exception
        raise AssertionError('wrong exception type')

with assert_raises(KeyError):
    {}['foo']

使用生成器定义上下文管理器

当讨论生成器时，据说我们相比实现为类的迭代器更倾向于生成器，因为它们更短小方便，状态被局部保存而非实例和变量中。另一方面，正如双向通信章节描述的那样，生成器和它的调用者之间的数据流可以是双向的。包括异常，可以直接传递给生成器。我们想将上下文管理器实现为特殊的生成器函数。事实上，生成器协议被设计成支持这个用例。

1
2
3
4
5
6
7

@contextlib.contextmanager
def some_generator(<arguments>):
    <setup>
    try:
        yield <value>
    finally:
        <cleanup>

contextlib.contextmanager装饰一个生成器并转换为上下文管理器。生成器必须遵循一些被包装(wrapper)函数强制执行的法则——最重要的是它至少yield一次。yield之前的部分从__enter__执行，上下文管理器中的代码块当生成器停在yield时执行，剩下的在__exit__中执行。如果异常被抛出，解释器通过__exit__的参数将之传递给包装函数，包装函数于是在yield语句处抛出异常。通过使用生成器，上下文管理器变得更短小精炼。

让我们用生成器重写closing的例子：

1
2
3
4
5
6

@contextlib.contextmanager
def closing(obj):
    try:
        yield obj
    finally:
        obj.close()

再把assert_raises改写成生成器：

1
2
3
4
5
6
7
8
9
10

@contextlib.contextmanager
def assert_raises(type):
    try:
        yield
    except type:
        return
    except Exception as value:
        raise AssertionError('wrong exception type')
    else:
        raise AssertionError('exception expected')

这里我们用装饰器将生成函数转化为上下文管理器！

FootNotes

[^2]:好吧它已经发布了= =，虽然在大多linux发行版中还是2.x和3.2。

[^5]:云里雾里= =

快速ebuild向导

这页是个非常简洁的ebuild写作指南。它不包含许多开发者面临的细节和问题，而是给出足够用来理解ebuild如何工作的微小的例子。

为了正确的涵盖所有来龙去脉，参见Ebuild Writing。General Concepts章节也很有用。

注意这里的例子，虽然基于真实的ebuilds树，但有些部分大刀阔斧的修剪、更改和简化了。

第一个Ebuild

我们从Exuberant Ctags工具开始，一个代码索引工具。这时一个简化的dev-util/ctags/ctags-5.5.4.ebuild(你能在主目录下找到真实的ebuild)

# Copyright 1999-2012 Gentoo Foundation
# Distributed under the terms of the GNU General Public License v2
# $Header: $

EAPI=4

DESCRIPTION="Exuberant ctags generates tags files for quick source navigation"
HOMEPAGE="http://ctags.sourceforge.net"
SRC_URI="mirror://sourceforge/ctags/${P}.tar.gz"

LICENSE="GPL-2"
SLOT="0"
KEYWORDS="~mips ~sparc ~x86"

src_configure() {
    econf --with-posix-regex
}

src_install() {
    emake DESTDIR="${D}" install

    dodoc FAQ NEWS README
    dohtml EXTENDING.html ctags.html
}

基本格式

如你所见，ebuilds仅仅是在特殊环境中执行的bash脚本。

在ebuild的顶部是头块(header block)，出现在所有ebuild中。

Ebuild使用tabs缩进，每个tab代表四格空格。参见Ebuild File Format.

信息变量

接着，有一系列变量将告诉Portage有关包和ebuild的各种东西。

ebuild的EAPI，参见EAPI Usage and Description

DESCRIPTION变量是包及包的作用的简短描述。

HOMEPAGE是链接到包的主页的链接。(切记包含http://部分)。

SRC_URI告诉Portage用来下载源码包的地址。这里，mirror://sourceforge/是意为“任何Sourceforge镜像”的特殊标记。${P}是由Portage设置的只读变量，即包名和版本——示例中是ctags-5.5.4。

LICENSE是协议GPL-2(GNU General Public License version 2)。

SLOT告诉Portage这个包安装到哪个slot。

KEYWORDS变量设置ebuild测试的架构。我们使用～keyword给我们新写的ebuild，包不允许被直接推送到稳定版，即使他们似乎工作。参见Keywording查看细节。

构建函数

接着一个叫src_configure的函数，Portage将在配置(configure)包时调用它。econf是执行./configure的一个封装。如果由于某些原因在econf时出错，Portage将停止而非继续安装。

当Portage准备安装包时会调用src_install函数。这里有点微妙——并非直接安装到文件系统，我们必须安装到一个Portage通过${D}变量给出的特殊位置(Portage设置这个——参见Install Destinations和Sandbox)

注意：常规安装方法是emake DESTDIR="${D}" install，这适合所有符合标准的Makefile。如果给出sandbox错误，尝试用einstall代替。如果仍然失败，参看src_install如何手工安装。

dodoc和dohtml部分是安装文件到相应的/usr/share/doc部分的辅助函数。

ebuild可以定义其它函数(参见Ebuild Functions)。在大多情况下，Portage提供合理的默认实现，通常做正确的事情，不需要定义src_unpack和src_compile函数。例如，src_unpack函数被用来解包或给源码打补丁，但是这个例子中默认实现做了我们所需要的所有事情。同样默认的src_compile函数将调用emake——一个make的封装。

注意：先前|| die结构不得不添加到每个命令后去检查错误。这在EAPI 4中不在必要——如果什么出错的话Portage提供的函数将自己die。

含依赖的ebuild

在ctags的例子中，我们没告诉Portage有关任何依赖。当情况是这样时，没关系，因为ctags仅仅需要一个基本的工具链来编译和运行(参见Implicit System Dependency理解为何我们不需要显式依赖)。然而事情很少这么简单。

这是app-misc/detox/detox-1.1.1.ebuild：

# Copyright 1999-2012 Gentoo Foundation
# Distributed under the terms of the GNU General Public License v2
# $Header: $

EAPI=4

DESCRIPTION="detox safely removes spaces and strange characters from filenames"
HOMEPAGE="http://detox.sourceforge.net/"
SRC_URI="mirror://sourceforge/${PN}/${P}.tar.bz2"

LICENSE="BSD"
SLOT="0"
KEYWORDS="~hppa ~mips sparc x86"

RDEPEND="dev-libs/popt"
DEPEND="${RDEPEND}
    sys-devel/flex
    sys-devel/bison"

src_configure() {
    econf --with-popt
}

src_install() {
    emake DESTDIR="${D}" install
    dodoc README CHANGES
}

你再次看到ebuild头和不通的信息变量。在SRC_URI中，${PN}用来获取不含尾部版本的包名(还有更多的这种变量——参见Predefined Read-Only Variables)。

我们再次定义了src_configure和src_install函数。

Portage依靠DEPEND和RDEPEND变量决定构建和运行包需要哪些变量。DEPEND变量列出编译时依赖，RDEPEND变量列出运行时依赖。参见Dependencies获取更复杂的例子。

带补丁的ebuild

我们经常要打补丁。这通过epatch辅助函数在src_prepare函数中完成。为了使用epatch，必须告诉Portage需要的eutilseclass(eclass就像库一样)——这通过在ebuild顶部的inherit eutils完成。这是app-misc/detoxdetox-1.1.0.ebuild：

# Copyright 1999-2012 Gentoo Foundation
# Distributed under the terms of the GNU General Public License v2
# $Header: $

EAPI=4

inherit eutils

DESCRIPTION="detox safely removes spaces and strange characters from filenames"
HOMEPAGE="http://detox.sourceforge.net/"
SRC_URI="mirror://sourceforge/${PN}/${P}.tar.bz2"

LICENSE="BSD"
SLOT="0"
KEYWORDS="~hppa ~mips ~sparc ~x86"

RDEPEND="dev-libs/popt"
DEPEND="${RDEPEND}
    sys-devel/flex
    sys-devel/bison"

src_prepare() {
    epatch "${FILESDIR}"/${P}-destdir.patch \
        "${FILESDIR}"/${P}-parallel_build.patch
}

src_configure() {
    econf --with-popt
}

src_install() {
    emake DESTDIR="${D}" install
    dodoc README CHANGES
}

注意${FILESDIR}/${P}-destdir.patch指向detox-1.1.0-destdir.patch，这个文件在Portage树的files/子文件夹中。更大的补丁文件必须在你的开发者空间dev.gentoo.org而不是files/或镜像中——参见Gentoo Mirrors和Patching with epatch。

带USE标记的ebuild

对USE标记，这里有个例子dev-libs/libiconv/libiconv-1.9.2.ebuild，一个libc实现的iconv替代。

# Copyright 1999-2012 Gentoo Foundation
# Distributed under the terms of the GNU General Public License v2
# $Header: $

EAPI=4

DESCRIPTION="GNU charset conversion library for libc which doesn't implement it"
HOMEPAGE="http://www.gnu.org/software/libiconv/"
SRC_URI="ftp://ftp.gnu.org/pub/gnu/libiconv/${P}.tar.gz"

LICENSE="LGPL-2.1"
SLOT="0"
KEYWORDS="~amd64 ~ppc ~sparc ~x86"
IUSE="nls"

DEPEND="!sys-libs/glibc"

src_configure() {
    econf $(use_enable nls)
}

src_install() {
    emake DESTDIR="${D}" install
}

注意IUSE变量。这列出了所有(非特殊)ebuild使用的use标记。除其它事项外，它还将被用作emerge -pv时输出。

这个包的./configure脚本使用了常规的--enable-nls或--disable-nls参数。我们用use_enable工具函数依赖用户USE标记自动生成这个(参见Query Function Reference)。

另一个复杂的例子是mail-client/sylpheed/sylpheed-1.0.4.ebuild：

# Copyright 1999-2012 Gentoo Foundation
# Distributed under the terms of the GNU General Public License v2
# $Header: $

EAPI=4

inherit eutils

DESCRIPTION="A lightweight email client and newsreader"
HOMEPAGE="http://sylpheed.good-day.net/"
SRC_URI="mirror://sourceforge/${PN}/${P}.tar.bz2"

LICENSE="GPL-2"
SLOT="0"
KEYWORDS="alpha amd64 hppa ia64 ppc ppc64 sparc x86"
IUSE="crypt imlib ipv6 ldap nls pda ssl xface"

RDEPEND="=x11-libs/gtk+-2*
    crypt? ( >=app-crypt/gpgme-0.4.5 )
    imlib? ( media-libs/imlib2 )
    ldap? ( >=net-nds/openldap-2.0.11 )
    pda? ( app-pda/jpilot )
    ssl? ( dev-libs/openssl )
    xface? ( >=media-libs/compface-1.4 )
    app-misc/mime-types
    x11-misc/shared-mime-info"
DEPEND="${RDEPEND}
    dev-util/pkgconfig
    nls? ( >=sys-devel/gettext-0.12.1 )"

src_prepare() {
    epatch "${FILESDIR}"/${PN}-namespace.diff \
        "${FILESDIR}"/${PN}-procmime.diff
}

src_configure() {
    econf \
        $(use_enable nls) \
        $(use_enable ssl) \
        $(use_enable crypt gpgme) \
        $(use_enable pda jpilot) \
        $(use_enable ldap) \
        $(use_enable ipv6) \
        $(use_enable imlib) \
        $(use_enable xface compface)
}

src_install() {
    emake DESTDIR="${D}" install

    doicon sylpheed.png
    domenu sylpheed.desktop

    dodoc [A-Z][A-Z]* ChangeLog*
}

注意可选依赖。有些use_enable行使用两个参数的版本——这在USE标记名不完全匹配./configure参数时非常有用。

Image to CSS

Warning： 我警告过了，浏览器弱爆了的就不要点开了= =。内存太小了就慎重打开= =

眼残，竟然看到这种奇葩的演示

看了看CSS源码，想起来可以用Python Image Library做一个generator。

结果就真得做了一个……

Codes

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

#! /bin/env python
# -*- coding: utf-8 -*-

"""
Script to turn image into css
"""

import Image
import sys

__author__ = "Reverland (lhtlyy@gmail.com)"


def getcss(im):
    """docstring for get"""
    css = """position: absolute;
    top: 30px;
    left: 30px;
    width: 0;
    height: 0;
    box-shadow:
        """
    string = '%dpx %dpx 0px 1px rgb%s,\n'
    for y in range(0, im.size[1], 1):
        for x in range(0, im.size[0], 1):
            if im.size[1] - y <= 1 and im.size[0] - x <= 1:
                string = '%dpx %dpx 0px 1px rgb%s;\n'
            color = im.getpixel((x, y))
            css += string % (x, y, color)
    return css


def gethtml(css):
    """docstring for gethtml"""
    html = """
    <div style="
    %s"></div>
    """ % css
    return html

if __name__ == '__main__':
    filename = sys.argv[1]
    #outfile = sys.argv[2]
    im = Image.open(filename)
    ratio = 0.5
    size = (int(ratio * im.size[0]), int(ratio * im.size[1]))
    im.thumbnail(size)
    html = gethtml(getcss(im))
    print html
    # with open(outfile, 'wb') as f:
    #     f.write(html)

Demo

点击显示图像

意义？

意义在于可能你会死机= =

尝试复制看看？

Python小练习：追踪导弹仿真

警告：浏览器可能不支持html嵌入标签，那么很抱歉什么视频也看不到。建议使用最新稳定版的firefox/chromium

仿真的时候面向对象会很方便。

缘起

某天终于没有雾霾的时候，我在操场上追赶正在散步的父亲……然后，我就想多了= =跟踪导弹是个什么轨迹呢？

仿真

首先把问题简化下，如果从圆心开始追赶圆上匀速运动的物体，是什么情况。

我先自己设法用微分笔算了算，发现实在搞不定。

上网查查导弹问题看到一些简单的直线问题，都涉及一堆微分方程和欧拉法迭代啥的……

干脆自己仿真下吧。这是最初的版本，完全没有面向对象概念。

前半部分调整图像的代码完全可以不看，从while循环开始即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

import matplotlib.pyplot as plt
import numpy as np
tolerance = 1e-1
radius = np.pi
v_o = 20
x_o, y_o = 0, radius

x_m, y_m = -radius, 0
v_m = 5

plt.figure(figsize=(10, 10), dpi=80)
plt.title(" missile flight simulator ", fontsize=40)
plt.xlim(-4, 4)
plt.ylim(-4, 4)
#plt.xticks([])
#plt.yticks([])

# set spines
ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data', 0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data', 0))
plt.xticks([-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi], [r'$-\pi$', r'$-\pi/2$', r'$0$', r'$+\pi/2$', r'$+\pi$'])
plt.yticks([-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi], [r'$-\pi$', r'$-\pi/2$', r'$0$', r'$+\pi/2$', r'$+\pi$'])

# Note object and missile
plt.annotate('object start point', xy=(x_o, y_o),  xycoords='data',
             xytext=(+15, +15), textcoords='offset points', fontsize=12,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate('missile start point', xy=(x_m, y_m),  xycoords='data',
             xytext=(+15, +15), textcoords='offset points', fontsize=12,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))

# alpha labels
for label in ax.get_xticklabels() + ax.get_yticklabels():
    label.set_fontsize(16)
    label.set_bbox(dict(facecolor='white', edgecolor='None', alpha=0.65))


while True:
    if x_o == 0 and y_o == radius:
        beta = 0
    elif x_o == 0 and y_o == radius:
        beta = np.pi
    elif x_o < 0:
        beta = np.pi / 2 * 3 - np.arctan(y_o / x_o)
    else:
        beta = np.pi / 2 - np.arctan(y_o / x_o)
    if np.sqrt((x_o - x_m) ** 2 + (y_o - y_m) ** 2) < tolerance:
        print "collision"
        plt.plot(x_m, y_m, 'o')
        plt.annotate('crash point', xy=(x_m, y_m),  xycoords='data',
                     xytext=(+15, +15), textcoords='offset points', fontsize=12,
                     arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
        plt.pause(0.1)
        break
    elif x_o < x_m:
        alpha = np.pi + np.arctan((y_o - y_m) / (x_o - x_m))
    elif x_o == x_m:
        alpha = np.pi / 2
    else:
        alpha = np.arctan((y_o - y_m) / (x_o - x_m))
    x_o = radius * np.sin(beta + v_o * 0.01 / np.pi / 2)
    y_o = radius * np.cos(beta + v_o * 0.01 / np.pi / 2)
    x_m = x_m + v_m * 0.01 * np.cos(alpha)
    y_m = y_m + v_m * 0.01 * np.sin(alpha)
    #print alpha, beta
    plt.plot(x_o, y_o, 'r.', alpha=.5)
    plt.plot(x_m, y_m, 'bx', alpha=.5)
    plt.legend(("target", "missile"), loc="upper left", prop={'size': 12})
    plt.pause(0.1)

我发现如果我速度够慢，未必追得上，甚至连被追踪物的轨道都不会进入……这挺出乎意料的，本来以为一定追的上

回到最初的问题，我从我的位置上追我父亲（好傻，都不知道估计下位置……）

面向对象

问题来了，如果我要仿真不只一个追踪导弹，比如还想仿真一个拦截导弹呢？

拦截失败演示

还可以用上面的方法不断扩充代码，每个对象写重复的代码。

但这时面向对象就能发挥威力，减少代码重用了。

以下是对四蜗牛聚合线问题的仿真。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

import matplotlib.pyplot as plt
import numpy as np
tolerance = 1e-1
radius = np.pi

# missile 1
x_m1, y_m1 = -np.pi, 0
v_m1 = 5

# missile 2
x_m2, y_m2 = 0, np.pi
v_m2 = v_m1
# missile 3
x_m3, y_m3 = np.pi, 0
v_m3 = v_m1
# missile 4
x_m4, y_m4 = 0, -np.pi
v_m4 = v_m1

plt.figure(figsize=(10, 10), dpi=80)
plt.title(" missile flight simulator ", fontsize=40)
plt.xlim(-4, 4)
plt.ylim(-4, 4)
#plt.xticks([])
#plt.yticks([])

# set spines
ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data', 0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data', 0))
plt.xticks([-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi], [r'$-\pi$', r'$-\pi/2$', r'$0$', r'$+\pi/2$', r'$+\pi$'])
plt.yticks([-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi], [r'$-\pi$', r'$-\pi/2$', r'$0$', r'$+\pi/2$', r'$+\pi$'])

plt.annotate('missile start point', xy=(x_m1, y_m1),  xycoords='data',
             xytext=(+15, +15), textcoords='offset points', fontsize=12,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))

# alpha labels
for label in ax.get_xticklabels() + ax.get_yticklabels():
    label.set_fontsize(16)
    label.set_bbox(dict(facecolor='white', edgecolor='None', alpha=0.65))


class ob(object):
    """docstring for ob"""
    def __init__(self, x, y):
        self.x = x
        self.y = y


class missile(ob):
    """docstring for missile"""
    def __init__(self, x, y):
        super(missile, self).__init__(x, y)

    def forward(self, v, target):
        """docstring for forward"""
        if self.x < target.x:
            alpha = np.arctan((target.y - self.y) / (target.x - self.x))
        elif self.x > target.x:
            alpha = np.pi + np.arctan((target.y - self.y) / (target.x - self.x))
        elif self.x == target.x and self.y < target.y:
            alpha = np.pi / 2
        else:
            alpha = -np.pi / 2
        self.x = self.x + v * 0.01 * np.cos(alpha)
        self.y = self.y + v * 0.01 * np.sin(alpha)
        return self.x, self.y

    def distance(self, target):
        """docstring for distance"""
        return np.sqrt((self.x - target.x) ** 2 + (self.y - target.y) ** 2)


class target(ob):
    """docstring for target"""
    def __init__(self, x, y):
        super(target, self).__init__(x, y)

    def newposition(self, x, y):
        """docstring for newposition"""
        self.x = x
        self.y = y

m1 = missile(x_m1, y_m1)
m2 = missile(x_m2, y_m2)
m3 = missile(x_m3, y_m3)
m4 = missile(x_m4, y_m4)

while True:
    if m1.distance(m2) < tolerance or m1.distance(m3) < tolerance or m1.distance(m4) < tolerance:
        print "collision"
        plt.plot(x_m1, y_m1, 'o')
        plt.annotate('crash point', xy=(x_m1, y_m1),  xycoords='data',
                     xytext=(+15, +15), textcoords='offset points', fontsize=12,
                     arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
        plt.pause(0.1)
        plt.show()
        break
    elif m3.distance(m2) < tolerance or m3.distance(m4) < tolerance:
        print "collision"
        plt.plot(x_m3, y_m3, 'o')
        plt.annotate('crash point', xy=(x_m3, y_m3),  xycoords='data',
                     xytext=(+15, +15), textcoords='offset points', fontsize=12,
                     arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
        plt.pause(0.1)
        plt.show
        break
    x_m1, y_m1 = m1.forward(v_m1, m2)
    x_m2, y_m2 = m2.forward(v_m2, m3)
    x_m3, y_m3 = m3.forward(v_m3, m4)
    x_m4, y_m4 = m4.forward(v_m4, m1)
    #print alpha, beta
    plt.plot(x_m1, y_m1, 'bx', alpha=.5)
    plt.plot(x_m2, y_m2, 'k*', alpha=.5)
    plt.plot(x_m3, y_m3, 'r.', alpha=.5)
    plt.plot(x_m4, y_m4, 'gp', alpha=.5)
    plt.legend(("missile1", "missile2", "missile3", "missile4"), loc="upper left", prop={'size': 12})
    plt.pause(0.1)

总结

面向对象方法对仿真问题非常合适，能有效简化代码，做到DRY(Don’t repeat yourself)。

搞着玩的，也许我该想想复试怎么办了……

Generate Ascii Image Like the Matrix

Inspired by this online image to text converter and this post on Oschina

script hosts at github. Poor hack, Any advice is welcomed.

Edit: Oooops, Demo sucks on my site, it’s perfect if you generate an alone html file.look here and here

Demo